{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
       "@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
       "@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
       "@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
       "\n",
       "    div.cell{\n",
       "        width: 900px;\n",
       "        margin-left: 0% !important;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.text_cell code {\n",
       "        background: transparent;\n",
       "        color: #000000;\n",
       "        font-weight: 600;\n",
       "        font-size: 11pt;\n",
       "        font-style: bold;\n",
       "        font-family:  'Source Code Pro', Consolas, monocco, monospace;\n",
       "   }\n",
       "    h1 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "\t}\n",
       "\t\n",
       "    div.input_area {\n",
       "        background: #F6F6F9;\n",
       "        border: 1px solid #586e75;\n",
       "    }\n",
       "\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 30pt;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h2 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "        text-align: left;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 200;\n",
       "        font-size: 16pt;\n",
       "        font-style: italic;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1.5em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h3 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 200;\n",
       "        font-size: 14pt;\n",
       "        line-height: 100%;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 2em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 14pt;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 200;\n",
       "        font-style: normal;\n",
       "        color: #1d3b84;\n",
       "        font-size: 16pt;\n",
       "        margin-bottom: 0em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Fira sans', verdana,arial,sans-serif;\n",
       "        line-height: 125%;\n",
       "        font-size: 115%;\n",
       "        text-align:justify;\n",
       "        text-justify:inter-word;\n",
       "    }\n",
       "    div.output_subarea.output_text.output_pyout {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_subarea.output_stream.output_stdout.output_text {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_wrapper{\n",
       "        margin-top:0.2em;\n",
       "        margin-bottom:0.2em;\n",
       "}\n",
       "\n",
       "    code{\n",
       "      font-size: 70%;\n",
       "    }\n",
       "    .rendered_html code{\n",
       "    background-color: transparent;\n",
       "    }\n",
       "    ul{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li li{\n",
       "        padding-left: 0.2em; \n",
       "        margin-bottom: 0.2em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    ol{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ol li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    a:link{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:visited{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:hover{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:focus{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:active{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    .rendered_html :link {\n",
       "       text-decoration: underline; \n",
       "    }\n",
       "    .rendered_html :hover {\n",
       "       text-decoration: none; \n",
       "    }\n",
       "    .rendered_html :visited {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :focus {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :active {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    } \n",
       "    hr {\n",
       "      color: #f3f3f3;\n",
       "      background-color: #f3f3f3;\n",
       "      height: 1px;\n",
       "    }\n",
       "    blockquote{\n",
       "      display:block;\n",
       "      background: #fcfcfc;\n",
       "      border-left: 5px solid #c76c0c;\n",
       "      font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "      width:680px;\n",
       "      padding: 10px 10px 10px 10px;\n",
       "      text-align:justify;\n",
       "      text-justify:inter-word;\n",
       "      }\n",
       "      blockquote p {\n",
       "        margin-bottom: 0;\n",
       "        line-height: 125%;\n",
       "        font-size: 100%;\n",
       "      }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    scale:100,\n",
       "                        availableFonts: [],\n",
       "                        preferredFont:null,\n",
       "                        webFont: \"TeX\",\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# format the book\n",
    "%matplotlib inline\n",
    "import sys\n",
    "from __future__ import division, print_function\n",
    "import sys\n",
    "sys.path.insert(0,'../code')\n",
    "import book_format\n",
    "book_format.load_style('../code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dealing with Dimensions\n",
    "\n",
    "## Belly button bacteria\n",
    "\n",
    "Belly Button Biodiversity 2.0 (BBB2) is a nation-wide citizen science\n",
    "project with the goal of identifying bacterial species that can be found\n",
    "in human navels (<http://bbdata.yourwildlife.org>). The project might\n",
    "seem whimsical, but it is part of an increasing interest in the human\n",
    "microbiome, the set of microorganisms that live on human skin and parts\n",
    "of the body.\n",
    "\n",
    "In their pilot study, BBB2 researchers collected swabs from the navels\n",
    "of 60 volunteers, used multiplex pyrosequencing to extract and sequence\n",
    "fragments of 16S rDNA, then identified the species or genus the\n",
    "fragments came from. Each identified fragment is called a “read.”\n",
    "\n",
    "We can use these data to answer several related questions:\n",
    "\n",
    "-   Based on the number of species observed, can we estimate the total\n",
    "    number of species in the environment?\n",
    "\n",
    "-   Can we estimate the prevalence of each species; that is, the\n",
    "    fraction of the total population belonging to each species?\n",
    "\n",
    "-   If we are planning to collect additional samples, can we predict how\n",
    "    many new species we are likely to discover?\n",
    "\n",
    "-   How many additional reads are needed to increase the fraction of\n",
    "    observed species to a given threshold?\n",
    "\n",
    "These questions make up what is called the **Unseen Species\n",
    "problem**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lions and tigers and bears\n",
    "\n",
    "I’ll start with a simplified version of the problem where we know that\n",
    "there are exactly three species. Let’s call them lions, tigers and\n",
    "bears. Suppose we visit a wild animal preserve and see 3 lions, 2 tigers\n",
    "and one bear.\n",
    "\n",
    "If we have an equal chance of observing any animal in the preserve, the\n",
    "number of each species we see is governed by the multinomial\n",
    "distribution. If the prevalence of lions and tigers and bears is\n",
    "`p_lion` and `p_tiger` and `p_bear`, the likelihood of seeing 3 lions, 2\n",
    "tigers and one bear is\n",
    "\n",
    "```python\n",
    "    p_lion**3 * p_tiger**2 * p_bear**1\n",
    "```\n",
    "\n",
    "An approach that is tempting, but not correct, is to use beta\n",
    "distributions, as in Section [beta], to describe the prevalence of each\n",
    "species separately. For example, we saw 3 lions and 3 non-lions; if we\n",
    "think of that as 3 “heads” and 3 “tails,” then the posterior\n",
    "distribution of `p_lion` is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'Beta' object has no attribute 'MaximumLikelihood'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-3-ca1777461973>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mbeta\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mthinkbayes\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mBeta\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mbeta\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mUpdate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m3\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbeta\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMaximumLikelihood\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m: 'Beta' object has no attribute 'MaximumLikelihood'"
     ]
    }
   ],
   "source": [
    "import thinkbayes\n",
    "\n",
    "beta = thinkbayes.Beta()\n",
    "beta.Update((3, 3))\n",
    "print(beta.MaximumLikelihood())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maximum likelihood estimate for `p_lion` is the observed rate, 50%.\n",
    "Similarly the MLEs for `p_tiger` and `p_bear` are 33% and 17%.\n",
    "\n",
    "But there are two problems:\n",
    "\n",
    "1.  We have implicitly used a prior for each species that is uniform\n",
    "    from 0 to 1, but since we know that there are three species, that\n",
    "    prior is not correct. The right prior should have a mean of 1/3, and\n",
    "    there should be zero likelihood that any species has a prevalence of\n",
    "    100%.\n",
    "\n",
    "2.  The distributions for each species are not independent, because the\n",
    "    prevalences have to add up to 1. To capture this dependence, we need\n",
    "    a joint distribution for the three prevalences.\n",
    "\n",
    "We can use a Dirichlet distribution to solve both of these problems (see\n",
    "<http://en.wikipedia.org/wiki/Dirichlet_distribution>). In the same way\n",
    "we used the beta distribution to describe the distribution of bias for a\n",
    "coin, we can use a Dirichlet distribution to describe the joint\n",
    "distribution of `p_lion`, `p_tiger` and `p_bear`.\n",
    "\n",
    "The Dirichlet distribution is the multi-dimensional generalization of\n",
    "the beta distribution. Instead of two possible outcomes, like heads and\n",
    "tails, the Dirichlet distribution handles any number of outcomes: in\n",
    "this example, three species.\n",
    "\n",
    "If there are `n` outcomes, the Dirichlet distribution is\n",
    "described by `n` parameters, written $\\alpha_1$ through\n",
    "$\\alpha_n$.\n",
    "\n",
    "Here’s the definition, from `thinkbayes.py`, of a class that\n",
    "represents a Dirichlet distribution:\n",
    "\n",
    "```python\n",
    "    class Dirichlet(object):\n",
    "\n",
    "        def __init__(self, n):\n",
    "            self.n = n\n",
    "            self.params = numpy.ones(n, dtype=numpy.int)\n",
    "```\n",
    "\n",
    "`n` is the number of dimensions; initially the parameters are\n",
    "all 1. I use a `numpy` array to store the parameters so I can\n",
    "take advantage of array operations.\n",
    "\n",
    "Given a Dirichlet distribution, the marginal distribution for each\n",
    "prevalence is a beta distribution, which we can compute like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MarginalBeta(self, i):\n",
    "    alpha0 = self.params.sum()\n",
    "    alpha = self.params[i]\n",
    "    return Beta(alpha, alpha0-alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`i` is the index of the marginal distribution we want.\n",
    "`alpha0` is the sum of the parameters; `alpha` is\n",
    "the parameter for the given species.\n",
    "\n",
    "In the example, the prior marginal distribution for each species is\n",
    "`Beta(1, 2)`. We can compute the prior means like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.333333333333\n",
      "0.333333333333\n",
      "0.333333333333\n"
     ]
    }
   ],
   "source": [
    "dirichlet = thinkbayes.Dirichlet(3)\n",
    "for i in range(3):\n",
    "    beta = dirichlet.MarginalBeta(i)\n",
    "    print(beta.Mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the prior mean prevalence for each species is 1/3.\n",
    "\n",
    "To update the Dirichlet distribution, we add the observations to the\n",
    "parameters like this:\n",
    "\n",
    "```python\n",
    "    def Update(self, data):\n",
    "        m = len(data)\n",
    "        self.params[:m] += data\n",
    "```\n",
    "\n",
    "`data` is a sequence of counts in the same order as\n",
    "`params`, so in this example, it should be the number of\n",
    "lions, tigers and bears.\n",
    "\n",
    "`data` can be shorter than `params`; in that case\n",
    "there are some species that have not been observed.\n",
    "\n",
    "Here’s code that updates `dirichlet` with the observed data\n",
    "and computes the posterior marginal distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.444444434084\n",
      "1 0.333333370657\n",
      "2 0.222325930766\n"
     ]
    }
   ],
   "source": [
    "data = [3, 2, 1]\n",
    "dirichlet.Update(data)\n",
    "\n",
    "for i in range(3):\n",
    "    beta = dirichlet.MarginalBeta(i)\n",
    "    pmf = beta.MakePmf()\n",
    "    print(i, pmf.Mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Distribution of prevalences for three species.](figs/species1.pdf)\n",
    "\n",
    "[fig.species1]\n",
    "\n",
    "Figure [fig.species1] shows the results. The posterior mean prevalences\n",
    "are 44%, 33%, and 22%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The hierarchical version\n",
    "\n",
    "We have solved a simplified version of the problem: if we know how many\n",
    "species there are, we can estimate the prevalence of each.\n",
    "\n",
    "Now let’s get back to the original problem, estimating the total number\n",
    "of species. To solve this problem I’ll define a meta-Suite, which is a\n",
    "Suite that contains other Suites as hypotheses. In this case, the\n",
    "top-level Suite contains hypotheses about the number of species; the\n",
    "bottom level contains hypotheses about prevalences.\n",
    "\n",
    "Here’s the class definition:\n",
    "\n",
    "    class Species(thinkbayes.Suite):\n",
    "\n",
    "        def __init__(self, ns):\n",
    "            hypos = [thinkbayes.Dirichlet(n) for n in ns]\n",
    "            thinkbayes.Suite.__init__(self, hypos)\n",
    "\n",
    "`__init__` takes a list of possible values for `n` and makes\n",
    "a list of Dirichlet objects.\n",
    "\n",
    "Here’s the code that creates the top-level suite:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from species import Species\n",
    "\n",
    "ns = range(3, 30)\n",
    "suite = Species(ns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`ns` is the list of possible values for `n`. We\n",
    "have seen 3 species, so there have to be at least that many. I chose an\n",
    "upper bound that seems reasonable, but we will check later that the\n",
    "probability of exceeding this bound is low. And at least initially we\n",
    "assume that any value in this range is equally likely.\n",
    "\n",
    "To update a hierarchical model, you have to update all levels. Usually\n",
    "you have to update the bottom level first and work up, but in this case\n",
    "we can update the top level first:\n",
    "\n",
    "```python\n",
    "#class Species\n",
    "\n",
    "    def Update(self, data):\n",
    "        thinkbayes.Suite.Update(self, data)\n",
    "        for hypo in self.Values():\n",
    "            hypo.Update(data)\n",
    "```\n",
    "\n",
    "`Species.Update` invokes `Update` in the parent\n",
    "class, then loops through the sub-hypotheses and updates them.\n",
    "\n",
    "Now all we need is a likelihood function:\n",
    "\n",
    "```python\n",
    "# class Species\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        dirichlet = hypo\n",
    "        like = 0\n",
    "        for i in range(1000):\n",
    "            like += dirichlet.Likelihood(data)\n",
    "\n",
    "        return like\n",
    "```\n",
    "\n",
    "`data` is a sequence of observed counts; `hypo` is\n",
    "a Dirichlet object. `Species.Likelihood` calls\n",
    "`Dirichlet.Likelihood` 1000 times and returns the total.\n",
    "\n",
    "Why call it 1000 times? Because `Dirichlet.Likelihood`\n",
    "doesn’t actually compute the likelihood of the data under the whole\n",
    "Dirichlet distribution. Instead, it draws one sample from the\n",
    "hypothetical distribution and computes the likelihood of the data under\n",
    "the sampled set of prevalences.\n",
    "\n",
    "Here’s what it looks like:\n",
    "\n",
    "```python\n",
    "# class Dirichlet\n",
    "\n",
    "    def Likelihood(self, data):\n",
    "        m = len(data)\n",
    "        if self.n < m:\n",
    "            return 0\n",
    "\n",
    "        x = data\n",
    "        p = self.Random()\n",
    "        q = p[:m]**x\n",
    "        return q.prod()\n",
    "```\n",
    "\n",
    "The length of `data` is the number of species observed. If we\n",
    "see more species than we thought existed, the likelihood is 0.\n",
    "\n",
    "Otherwise we select a random set of prevalences, `p`, and\n",
    "compute the multinomial PMF, which is\n",
    "\n",
    "$$c_x  p_1^{x_1} \\cdots p_n^{x_n}$$ \n",
    "\n",
    "$p_i$ is the prevalence of the $i$th\n",
    "species, and $x_i$ is the observed number. The first term, $c_x$, is the\n",
    "multinomial coefficient; I leave it out of the computation because it is\n",
    "a multiplicative factor that depends only on the data, not the\n",
    "hypothesis, so it gets normalized away (see\n",
    "<http://en.wikipedia.org/wiki/Multinomial_distribution>).\n",
    "\n",
    "`m` is the number of observed species. We only need the first\n",
    "`m` elements of `p`; for the others, $x_i$ is 0,\n",
    "so $p_i^{x_i}$ is 1, and we can leave them out of the product."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random sampling\n",
    "\n",
    "There are two ways to generate a random sample from a Dirichlet\n",
    "distribution. One is to use the marginal beta distributions, but in that\n",
    "case you have to select one at a time and scale the rest so they add up\n",
    "to 1 (see\n",
    "<http://en.wikipedia.org/wiki/Dirichlet_distribution#Random_number_generation>).\n",
    "\n",
    "A less obvious, but faster, way is to select values from `n`\n",
    "gamma distributions, then normalize by dividing through by the total.\n",
    "Here’s the code:\n",
    "\n",
    "```python\n",
    "# class Dirichlet\n",
    "\n",
    "    def Random(self):\n",
    "        p = numpy.random.gamma(self.params)\n",
    "        return p / p.sum()\n",
    "```\n",
    "\n",
    "Now we’re ready to look at some results. Here is the code that extracts\n",
    "the posterior distribution of `n`:\n",
    "\n",
    " ```python\n",
    "  \n",
    "   def DistOfN(self):\n",
    "        pmf = thinkbayes.Pmf()\n",
    "        for hypo, prob in self.Items():\n",
    "            pmf.Set(hypo.n, prob)\n",
    "        return pmf\n",
    "```\n",
    "\n",
    "`DistOfN` iterates through the top-level hypotheses and\n",
    "accumulates the probability of each `n`.\n",
    "\n",
    "![Posterior distribution of `n`.](figs/species2.pdf)\n",
    "\n",
    "[fig.species2]\n",
    "\n",
    "Figure [fig.species2] shows the result. The most likely value is 4.\n",
    "Values from 3 to 7 are reasonably likely; after that the probabilities\n",
    "drop off quickly. The probability that there are 29 species is low\n",
    "enough to be negligible; if we chose a higher bound, we would get nearly\n",
    "the same result.\n",
    "\n",
    "Remember that this result is based on a uniform prior for\n",
    "`n`. If we have background information about the number of\n",
    "species in the environment, we might choose a different prior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization\n",
    "\n",
    "I have to admit that I am proud of this example. The Unseen Species\n",
    "problem is not easy, and I think this solution is simple and clear, and\n",
    "takes surprisingly few lines of code (about 50 so far).\n",
    "\n",
    "The only problem is that it is slow. It’s good enough for the example\n",
    "with only 3 observed species, but not good enough for the belly button\n",
    "data, with more than 100 species in some samples.\n",
    "\n",
    "The next few sections present a series of optimizations we need to make\n",
    "this solution scale. Before we get into the details, here’s a road map.\n",
    "\n",
    "-   The first step is to recognize that if we update the Dirichlet\n",
    "    distributions with the same data, the first `m`\n",
    "    parameters are the same for all of them. The only difference is the\n",
    "    number of hypothetical unseen species. So we don’t really need\n",
    "    `n` Dirichlet objects; we can store the parameters in the\n",
    "    top level of the hierarchy. `Species2` implements this\n",
    "    optimization.\n",
    "\n",
    "-   `Species2` also uses the same set of random values for\n",
    "    all of the hypotheses. This saves time generating random values, but\n",
    "    it has a second benefit that turns out to be more important: by\n",
    "    giving all hypotheses the same selection from the sample space, we\n",
    "    make the comparison between the hypotheses more fair, so it takes\n",
    "    fewer iterations to converge.\n",
    "\n",
    "-   Even with these changes there is a major performance problem. As the\n",
    "    number of observed species increases, the array of random\n",
    "    prevalences gets bigger, and the chance of choosing one that is\n",
    "    approximately right becomes small. So the vast majority of\n",
    "    iterations yield small likelihoods that don’t contribute much to the\n",
    "    total, and don’t discriminate between hypotheses.\n",
    "\n",
    "    The solution is to do the updates one species at a time.\n",
    "    `Species4` is a simple implementation of this strategy\n",
    "    using Dirichlet objects to represent the sub-hypotheses.\n",
    "\n",
    "-   Finally, `Species5` combines the sub-hypotheses into the\n",
    "    top level and uses `numpy` array operations to speed\n",
    "    things up.\n",
    "\n",
    "If you are not interested in the details, feel free to skip to\n",
    "Section [belly] where we look at results from the belly button data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Collapsing the hierarchy\n",
    "\n",
    "All of the bottom-level Dirichlet distributions are updated with the\n",
    "same data, so the first `m` parameters are the same for all\n",
    "of them. We can eliminate them and merge the parameters into the\n",
    "top-level suite. `Species2` implements this optimization:\n",
    "\n",
    "```python\n",
    "class Species2(object):\n",
    "\n",
    "    def __init__(self, ns):\n",
    "        self.ns = ns\n",
    "        self.probs = numpy.ones(len(ns), dtype=numpy.double)\n",
    "        self.params = numpy.ones(self.high, dtype=numpy.int)\n",
    "```\n",
    "\n",
    "`ns` is the list of hypothetical values for `n`;\n",
    "`probs` is the list of corresponding probabilities. And\n",
    "`params` is the sequence of Dirichlet parameters, initially\n",
    "all 1.\n",
    "\n",
    "`Species2.Update` updates both levels of the hierarchy: first\n",
    "the probability for each value of `n`, then the Dirichlet\n",
    "parameters:\n",
    "\n",
    "```python\n",
    "# class Species2\n",
    "\n",
    "    def Update(self, data):\n",
    "        like = numpy.zeros(len(self.ns), dtype=numpy.double)\n",
    "        for i in range(1000):\n",
    "            like += self.SampleLikelihood(data)\n",
    "\n",
    "        self.probs *= like\n",
    "        self.probs /= self.probs.sum()\n",
    "\n",
    "        m = len(data)\n",
    "        self.params[:m] += data\n",
    "```\n",
    "\n",
    "`SampleLikelihood` returns an array of likelihoods, one for\n",
    "each value of `n`. `like` accumulates the total\n",
    "likelihood for 1000 samples. `self.probs` is multiplied by\n",
    "the total likelihood, then normalized. The last two lines, which update\n",
    "the parameters, are the same as in `Dirichlet.Update`.\n",
    "\n",
    "Now let’s look at `SampleLikelihood`. There are two\n",
    "opportunities for optimization here:\n",
    "\n",
    "-   When the hypothetical number of species, `n`, exceeds the\n",
    "    observed number, `m`, we only need the first\n",
    "    `m` terms of the multinomial PMF; the rest are 1.\n",
    "\n",
    "-   If the number of species is large, the likelihood of the data might\n",
    "    be too small for floating-point (see  [underflow]). So it is safer\n",
    "    to compute log-likelihoods.\n",
    "\n",
    "Again, the multinomial PMF is \n",
    "\n",
    "$$c_x p_1^{x_1} \\cdots p_n^{x_n}$$ \n",
    "\n",
    "So the log-likelihood is\n",
    "\n",
    "$$\\log c_x + x_1 \\log p_1 + \\cdots + x_n \\log p_n$$\n",
    "\n",
    "which is fast and easy to compute. Again, $c_x$ it is the same for all\n",
    "hypotheses, so we can drop it. Here’s the code:\n",
    "\n",
    "```python\n",
    "# class Species2\n",
    "\n",
    "    def SampleLikelihood(self, data):\n",
    "        gammas = numpy.random.gamma(self.params)\n",
    "\n",
    "        m = len(data)\n",
    "        row = gammas[:m]\n",
    "        col = numpy.cumsum(gammas)\n",
    "\n",
    "        log_likes = []\n",
    "        for n in self.ns:\n",
    "            ps = row / col[n-1]\n",
    "            terms = data * numpy.log(ps)\n",
    "            log_like = terms.sum()\n",
    "            log_likes.append(log_like)\n",
    "\n",
    "        log_likes -= numpy.max(log_likes)\n",
    "        likes = numpy.exp(log_likes)\n",
    "\n",
    "        coefs = [thinkbayes.BinomialCoef(n, m) for n in self.ns]\n",
    "        likes *= coefs\n",
    "\n",
    "        return likes\n",
    "```\n",
    "\n",
    "`gammas` is an array of values from a gamma distribution; its\n",
    "length is the largest hypothetical value of `n`.\n",
    "`row` is just the first `m` elements of\n",
    "`gammas`; since these are the only elements that depend on\n",
    "the data, they are the only ones we need.\n",
    "\n",
    "For each value of `n` we need to divide `row` by\n",
    "the total of the first `n` values from `gamma`.\n",
    "`cumsum` computes these cumulative sums and stores them in\n",
    "`col`.\n",
    "\n",
    "The loop iterates through the values of `n` and accumulates a\n",
    "list of log-likelihoods.\n",
    "\n",
    "Inside the loop, `ps` contains the row of probabilities,\n",
    "normalized with the appropriate cumulative sum. `terms`\n",
    "contains the terms of the summation, $x_i \\log p_i$, and `log_like`\n",
    "contains their sum.\n",
    "\n",
    "After the loop, we want to convert the log-likelihoods to linear\n",
    "likelihoods, but first it’s a good idea to shift them so the largest\n",
    "log-likelihood is 0; that way the linear likelihoods are not too small\n",
    "(see  [underflow]).\n",
    "\n",
    "Finally, before we return the likelihood, we have to apply a correction\n",
    "factor, which is the number of ways we could have observed these\n",
    "`m` species, if the total number of species is\n",
    "`n`. `BinomialCoefficient` computes “n choose m”,\n",
    "which is written $\\binom{n}{m}$.\n",
    "\n",
    "As often happens, the optimized version is less readable and more\n",
    "error-prone than the original. But that’s one reason I think it is a\n",
    "good idea to start with the simple version; we can use it for regression\n",
    "testing. I plotted results from both versions and confirmed that they\n",
    "are approximately equal, and that they converge as the number of\n",
    "iterations increases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One more problem\n",
    "\n",
    "There’s more we could do to optimize this code, but there’s another\n",
    "problem we need to fix first. As the number of observed species\n",
    "increases, this version gets noisier and takes more iterations to\n",
    "converge on a good answer.\n",
    "\n",
    "The problem is that if the prevalences we choose from the Dirichlet\n",
    "distribution, the `ps`, are not at least approximately right,\n",
    "the likelihood of the observed data is close to zero and almost equally\n",
    "bad for all values of `n`. So most iterations don’t provide\n",
    "any useful contribution to the total likelihood. And as the number of\n",
    "observed species, `m`, gets large, the probability of\n",
    "choosing `ps` with non-negligible likelihood gets small.\n",
    "Really small.\n",
    "\n",
    "Fortunately, there is a solution. Remember that if you observe a set of\n",
    "data, you can update the prior distribution with the entire dataset, or\n",
    "you can break it up into a series of updates with subsets of the data,\n",
    "and the result is the same either way.\n",
    "\n",
    "For this example, the key is to perform the updates one species at a\n",
    "time. That way when we generate a random set of `ps`, only\n",
    "one of them affects the computed likelihood, so the chance of choosing a\n",
    "good one is much better.\n",
    "\n",
    "Here’s a new version that updates one species at a time:\n",
    "\n",
    "```python\n",
    "class Species4(Species):\n",
    "\n",
    "    def Update(self, data):\n",
    "        m = len(data)\n",
    "\n",
    "        for i in range(m):\n",
    "            one = numpy.zeros(i+1)\n",
    "            one[i] = data[i]            \n",
    "            Species.Update(self, one)\n",
    "```\n",
    "\n",
    "This version inherits `__init__` from `Species`, so it\n",
    "represents the hypotheses as a list of Dirichlet objects (unlike\n",
    "`Species2`).\n",
    "\n",
    "`Update` loops through the observed species and makes an\n",
    "array, `one`, with all zeros and one species count. Then it\n",
    "calls `Update` in the parent class, which computes the\n",
    "likelihoods and updates the sub-hypotheses.\n",
    "\n",
    "So in the running example, we do three updates. The first is something\n",
    "like “I have seen three lions.” The second is “I have seen two tigers\n",
    "and no additional lions.” And the third is “I have seen one bear and no\n",
    "more lions and tigers.”\n",
    "\n",
    "Here’s the new version of `Likelihood`:\n",
    "\n",
    "```python\n",
    "\n",
    "# class Species4\n",
    "\n",
    "    def Likelihood(self, data, hypo):\n",
    "        dirichlet = hypo\n",
    "        like = 0\n",
    "        for i in range(self.iterations):\n",
    "            like += dirichlet.Likelihood(data)\n",
    "\n",
    "        # correct for the number of unseen species the new one\n",
    "        # could have been\n",
    "        m = len(data)\n",
    "        num_unseen = dirichlet.n - m + 1\n",
    "        like *= num_unseen\n",
    "\n",
    "        return like\n",
    "```\n",
    "\n",
    "This is almost the same as `Species.Likelihood`. The\n",
    "difference is the factor, `num_unseen`. This correction is necessary\n",
    "because each time we see a species for the first time, we have to\n",
    "consider that there were some number of other unseen species that we\n",
    "might have seen. For larger values of `n` there are more\n",
    "unseen species that we could have seen, which increases the likelihood\n",
    "of the data.\n",
    "\n",
    "This is a subtle point and I have to admit that I did not get it right\n",
    "the first time. But again I was able to validate this version by\n",
    "comparing it to the previous versions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We’re not done yet\n",
    "\n",
    "Performing the updates one species at a time solves one problem, but it\n",
    "creates another. Each update takes time proportional to $k m$, where $k$\n",
    "is the number of hypotheses and $m$ is the number of observed species.\n",
    "So if we do $m$ updates, the total run time is proportional to $k m^2$.\n",
    "\n",
    "But we can speed things up using the same trick we used in\n",
    "Section [collapsing]: we’ll get rid of the Dirichlet objects and\n",
    "collapse the two levels of the hierarchy into a single object. So here’s\n",
    "yet another version of `Species`:\n",
    "\n",
    "```python\n",
    "class Species5(Species2):\n",
    "\n",
    "    def Update(self, data):\n",
    "        m = len(data)\n",
    "        for i in range(m):\n",
    "            self.UpdateOne(i+1, data[i])\n",
    "            self.params[i] += data[i]\n",
    "```\n",
    "\n",
    "This version inherits `__init__` from `Species2`, so it uses\n",
    "`ns` and `probs` to represent the distribution of\n",
    "`n`, and `params` to represent the parameters of\n",
    "the Dirichlet distribution.\n",
    "\n",
    "`Update` is similar to what we saw in the previous section.\n",
    "It loops through the observed species and calls `UpdateOne`:\n",
    "\n",
    "```python\n",
    "# class Species5\n",
    "\n",
    "    def UpdateOne(self, i, count):\n",
    "        likes = numpy.zeros(len(self.ns), dtype=numpy.double)\n",
    "        for i in range(self.iterations):\n",
    "            likes += self.SampleLikelihood(i, count)\n",
    "\n",
    "        unseen_species = [n-i+1 for n in self.ns]\n",
    "        likes *= unseen_species\n",
    "\n",
    "        self.probs *= likes\n",
    "        self.probs /= self.probs.sum()\n",
    "```\n",
    "\n",
    "This function is similar to `Species2.Update`, with two\n",
    "changes:\n",
    "\n",
    "-   The interface is different. Instead of the whole dataset, we get\n",
    "    `i`, the index of the observed species, and\n",
    "    `count`, how many of that species we’ve seen.\n",
    "\n",
    "-   We have to apply a correction factor for the number of unseen\n",
    "    species, as in `Species4.Likelihood`. The difference here\n",
    "    is that we update all of the likelihoods at once with array\n",
    "    multiplication.\n",
    "\n",
    "Finally, here’s `SampleLikelihood`:\n",
    "\n",
    "```python\n",
    "# class Species5\n",
    "\n",
    "    def SampleLikelihood(self, i, count):\n",
    "        gammas = numpy.random.gamma(self.params)\n",
    "\n",
    "        sums = numpy.cumsum(gammas)[self.ns[0]-1:]\n",
    "\n",
    "        ps = gammas[i-1] / sums\n",
    "        log_likes = numpy.log(ps) * count\n",
    "\n",
    "        log_likes -= numpy.max(log_likes)\n",
    "        likes = numpy.exp(log_likes)\n",
    "\n",
    "        return likes\n",
    "```\n",
    "\n",
    "This is similar to `Species2.SampleLikelihood`; the\n",
    "difference is that each update only includes a single species, so we\n",
    "don’t need a loop.\n",
    "\n",
    "The runtime of this function is proportional to the number of\n",
    "hypotheses, $k$. It runs $m$ times, so the run time of the update is\n",
    "proportional to $k m$. And the number of iterations we need to get an\n",
    "accurate result is usually small."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The belly button data\n",
    "\n",
    "That’s enough about lions and tigers and bears. Let’s get back to belly\n",
    "buttons. To get a sense of what the data look like, consider subject\n",
    "B1242, whose sample of 400 reads yielded 61 species with the following\n",
    "counts:\n",
    "\n",
    "    92, 53, 47, 38, 15, 14, 12, 10, 8, 7, 7, 5, 5, \n",
    "    4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2,\n",
    "    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
    "    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1\n",
    "\n",
    "There are a few dominant species that make up a large fraction of the\n",
    "whole, but many species that yielded only a single read. The number of\n",
    "these “singletons” suggests that there are likely to be at least a few\n",
    "unseen species.\n",
    "\n",
    "In the example with lions and tigers, we assume that each animal in the\n",
    "preserve is equally likely to be observed. Similarly, for the belly\n",
    "button data, we assume that each bacterium is equally likely to yield a\n",
    "read.\n",
    "\n",
    "In reality, each step in the data-collection process might introduce\n",
    "biases. Some species might be more likely to be picked up by a swab, or\n",
    "to yield identifiable amplicons. So when we talk about the prevalence of\n",
    "each species, we should remember this source of error.\n",
    "\n",
    "I should also acknowledge that I am using the term “species” loosely.\n",
    "First, bacterial species are not well defined. Second, some reads\n",
    "identify a particular species, others only identify a genus. To be more\n",
    "precise, I should say “operational taxonomic unit”, or OTU.\n",
    "\n",
    "Now let’s process some of the belly button data. I define a class called\n",
    "`Subject` to represent information about each subject in the\n",
    "study:\n",
    "\n",
    "```python\n",
    "class Subject(object):\n",
    "\n",
    "    def __init__(self, code):\n",
    "        self.code = code\n",
    "        self.species = []\n",
    "```\n",
    "\n",
    "Each subject has a string code, like “B1242”, and a list of (count,\n",
    "species name) pairs, sorted in increasing order by count.\n",
    "`Subject` provides several methods to make it easy to access\n",
    "these counts and species names. You can see the details in\n",
    "<http://thinkbayes.com/species.py>. For more information see\n",
    "Section [download].\n",
    "\n",
    "![Distribution of `n` for subject\n",
    "B1242.](figs/species-ndist-B1242.pdf)\n",
    "\n",
    "[species-ndist]\n",
    "\n",
    "`Subject` provides a method named `Process` that\n",
    "creates and updates a `Species5` suite, which represents the\n",
    "distributions of `n` and the prevalences.\n",
    "\n",
    "And `Suite2` provides `DistOfN`, which returns the\n",
    "posterior distribution of `n`.\n",
    "\n",
    "```python\n",
    "# class Suite2\n",
    "\n",
    "    def DistN(self):\n",
    "        items = zip(self.ns, self.probs)\n",
    "        pmf = thinkbayes.MakePmfFromItems(items)\n",
    "        return pmf\n",
    "```\n",
    "\n",
    "Figure [species-ndist] shows the distribution of `n` for\n",
    "subject B1242. The probability that there are exactly 61 species, and no\n",
    "unseen species, is nearly zero. The most likely value is 72, with 90%\n",
    "credible interval 66 to 79. At the high end, it is unlikely that there\n",
    "are as many as 87 species.\n",
    "\n",
    "Next we compute the posterior distribution of prevalence for each\n",
    "species. `Species2` provides `DistOfPrevalence`:\n",
    "\n",
    "    # class Species2\n",
    "\n",
    "        def DistOfPrevalence(self, index):\n",
    "            metapmf = thinkbayes.Pmf()\n",
    "\n",
    "            for n, prob in zip(self.ns, self.probs):\n",
    "                beta = self.MarginalBeta(n, index)\n",
    "                pmf = beta.MakePmf()\n",
    "                metapmf.Set(pmf, prob)\n",
    "\n",
    "            mix = thinkbayes.MakeMixture(metapmf)\n",
    "            return metapmf, mix\n",
    "\n",
    "`index` indicates which species we want. For each\n",
    "`n`, we have a different posterior distribution of\n",
    "prevalence.\n",
    "\n",
    "![Distribution of prevalences for subject\n",
    "B1242.](figs/species-prev-B1242.pdf)\n",
    "\n",
    "[species-prev]\n",
    "\n",
    "The loop iterates through the possible values of `n` and\n",
    "their probabilities. For each value of `n` it gets a Beta\n",
    "object representing the marginal distribution for the indicated species.\n",
    "Remember that Beta objects contain the parameters `alpha` and\n",
    "`beta`; they don’t have values and probabilities like a Pmf,\n",
    "but they provide `MakePmf`, which generates a discrete\n",
    "approximation to the continuous beta distribution.\n",
    "\n",
    "`metapmf` is a meta-Pmf that contains the distributions of\n",
    "prevalence, conditioned on `n`. `MakeMixture`\n",
    "combines the meta-Pmf into `mix`, which combines the\n",
    "conditional distributions into a single distribution of prevalence.\n",
    "\n",
    "Figure [species-prev] shows results for the five species with the most\n",
    "reads. The most prevalent species accounts for 23% of the 400 reads, but\n",
    "since there are almost certainly unseen species, the most likely\n",
    "estimate for its prevalence is 20%, with 90% credible interval between\n",
    "17% and 23%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predictive distributions\n",
    "\n",
    "![Simulated rarefaction curves for subject\n",
    "B1242.](figs/species-rare-B1242.pdf)\n",
    "\n",
    "[species-rare]\n",
    "\n",
    "I introduced the hidden species problem in the form of four related\n",
    "questions. We have answered the first two by computing the posterior\n",
    "distribution for `n` and the prevalence of each species.\n",
    "\n",
    "The other two questions are:\n",
    "\n",
    "-   If we are planning to collect additional reads, can we predict how\n",
    "    many new species we are likely to discover?\n",
    "\n",
    "-   How many additional reads are needed to increase the fraction of\n",
    "    observed species to a given threshold?\n",
    "\n",
    "To answer predictive questions like this we can use the posterior\n",
    "distributions to simulate possible future events and compute predictive\n",
    "distributions for the number of species, and fraction of the total, we\n",
    "are likely to see.\n",
    "\n",
    "The kernel of these simulations looks like this:\n",
    "\n",
    "1.  Choose `n` from its posterior distribution.\n",
    "\n",
    "2.  Choose a prevalence for each species, including possible unseen\n",
    "    species, using the Dirichlet distribution.\n",
    "\n",
    "3.  Generate a random sequence of future observations.\n",
    "\n",
    "4.  Compute the number of new species, `num_new`, as a function of the\n",
    "    number of additional reads, `k`.\n",
    "\n",
    "5.  Repeat the previous steps and accumulate the joint distribution of\n",
    "    `num_new` and `k`.\n",
    "\n",
    "And here’s the code. `RunSimulation` runs a single\n",
    "simulation:\n",
    "\n",
    "```python\n",
    "# class Subject\n",
    "\n",
    "    def RunSimulation(self, num_reads):\n",
    "        m, seen = self.GetSeenSpecies()\n",
    "        n, observations = self.GenerateObservations(num_reads)\n",
    "\n",
    "        curve = []\n",
    "        for k, obs in enumerate(observations):\n",
    "            seen.add(obs)\n",
    "\n",
    "            num_new = len(seen) - m\n",
    "            curve.append((k+1, num_new))\n",
    "\n",
    "        return curve\n",
    "```\n",
    "\n",
    "`num_reads` is the number of additional reads to simulate.\n",
    "`m` is the number of seen species, and `seen` is a\n",
    "set of strings with a unique name for each species. `n` is a\n",
    "random value from the posterior distribution, and\n",
    "`observations` is a random sequence of species names.\n",
    "\n",
    "Each time through the loop, we add the new observation to\n",
    "`seen` and record the number of reads and the number of new\n",
    "species so far.\n",
    "\n",
    "The result of `RunSimulation` is a **rarefaction\n",
    "curve**, represented as a list of pairs with the number of reads\n",
    "and the number of new species.\n",
    "\n",
    "Before we see the results, let’s look at `GetSeenSpecies` and\n",
    "`GenerateObservations`.\n",
    "\n",
    "```python\n",
    "#class Subject\n",
    "\n",
    "    def GetSeenSpecies(self):\n",
    "        names = self.GetNames()\n",
    "        m = len(names)\n",
    "        seen = set(SpeciesGenerator(names, m))\n",
    "        return m, seen\n",
    "```\n",
    "\n",
    "`GetNames` returns the list of species names that appear in\n",
    "the data files, but for many subjects these names are not unique. So I\n",
    "use `SpeciesGenerator` to extend each name with a serial\n",
    "number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def SpeciesGenerator(names, num):\n",
    "    i = 0\n",
    "    for name in names:\n",
    "        yield '%s-%d' % (name, i)\n",
    "        i += 1\n",
    "\n",
    "    while i < num:\n",
    "        yield 'unseen-%d' % i\n",
    "        i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a name like `Corynebacterium`,\n",
    "`SpeciesGenerator` yields `Corynebacterium-1`.\n",
    "When the list of names is exhausted, it yields names like\n",
    "`unseen-62`.\n",
    "\n",
    "Here is `GenerateObservations`:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a name like `Corynebacterium`,\n",
    "`SpeciesGenerator` yields `Corynebacterium-1`.\n",
    "When the list of names is exhausted, it yields names like\n",
    "`unseen-62`.\n",
    "\n",
    "Here is `GenerateObservations`:\n",
    "\n",
    "```python\n",
    "# class Subject\n",
    "\n",
    "    def GenerateObservations(self, num_reads):\n",
    "        n, prevalences = self.suite.SamplePosterior()\n",
    "\n",
    "        names = self.GetNames()\n",
    "        name_iter = SpeciesGenerator(names, n)\n",
    "\n",
    "        d = dict(zip(name_iter, prevalences))\n",
    "        cdf = thinkbayes.MakeCdfFromDict(d)\n",
    "        observations = cdf.Sample(num_reads)\n",
    "\n",
    "        return n, observations\n",
    "```\n",
    "\n",
    "Again, `num_reads` is the number of additional reads to generate.\n",
    "`n` and `prevalences` are samples from the\n",
    "posterior distribution.\n",
    "\n",
    "`cdf` is a Cdf object that maps species names, including the\n",
    "unseen, to cumulative probabilities. Using a Cdf makes it efficient to\n",
    "generate a random sequence of species names.\n",
    "\n",
    "Finally, here is `Species2.SamplePosterior`:\n",
    "\n",
    "```python\n",
    "    def SamplePosterior(self):\n",
    "        pmf = self.DistOfN()\n",
    "        n = pmf.Random()\n",
    "        prevalences = self.SamplePrevalences(n)\n",
    "        return n, prevalences\n",
    "```\n",
    "\n",
    "And `SamplePrevalences`, which generates a sample of\n",
    "prevalences conditioned on `n`:\n",
    "\n",
    "```python\n",
    "# class Species2\n",
    "\n",
    "    def SamplePrevalences(self, n):\n",
    "        params = self.params[:n]\n",
    "        gammas = numpy.random.gamma(params)\n",
    "        gammas /= gammas.sum()\n",
    "        return gammas\n",
    "```\n",
    "\n",
    "We saw this algorithm for generating random values from a Dirichlet\n",
    "distribution in Section [randomdir].\n",
    "\n",
    "Figure [species-rare] shows 100 simulated rarefaction curves for subject\n",
    "B1242. The curves are “jittered;” that is, I shifted each curve by a\n",
    "random offset so they would not all overlap. By inspection we can\n",
    "estimate that after 400 more reads we are likely to find 2–6 new\n",
    "species."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joint posterior\n",
    "\n",
    "![Distributions of the number of new species conditioned on the number\n",
    "of additional reads.](figs/species-cond-B1242.pdf)\n",
    "\n",
    "[species-cond]\n",
    "\n",
    "We can use these simulations to estimate the joint distribution of\n",
    "`num_new` and `k`, and from that we can get the distribution\n",
    "of `num_new` conditioned on any value of `k`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeJointPredictive(curves):\n",
    "    joint = thinkbayes.Joint()\n",
    "    for curve in curves:\n",
    "        for k, num_new in curve:\n",
    "            joint.Incr((k, num_new))\n",
    "    joint.Normalize()\n",
    "    return joint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`MakeJointPredictive` makes a Joint object, which is a Pmf\n",
    "whose values are tuples.\n",
    "\n",
    "`curves` is a list of rarefaction curves created by\n",
    "`RunSimulation`. Each curve contains a list of pairs of\n",
    "`k` and `num_new`.\n",
    "\n",
    "The resulting joint distribution is a map from each pair to its\n",
    "probability of occurring. Given the joint distribution, we can use\n",
    "`Joint.Conditional` get the distribution of `num_new`\n",
    "conditioned on `k` (see Section [conditional]).\n",
    "\n",
    "`Subject.MakeConditionals` takes a list of `ks`\n",
    "and computes the conditional distribution of `num_new` for each\n",
    "`k`. The result is a list of Cdf objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeConditionals(curves, ks):\n",
    "    joint = MakeJointPredictive(curves)\n",
    "\n",
    "    cdfs = []\n",
    "    for k in ks:\n",
    "        pmf = joint.Conditional(1, 0, k)\n",
    "        pmf.name = 'k=%d' % k\n",
    "        cdf = pmf.MakeCdf()\n",
    "        cdfs.append(cdf)\n",
    "\n",
    "    return cdfs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure [species-cond] shows the results. After 100 reads, the median\n",
    "predicted number of new species is 2; the 90% credible interval is 0 to\n",
    "5. After 800 reads, we expect to see 3 to 12 new species."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coverage\n",
    "\n",
    "![Complementary CDF of coverage for a range of additional\n",
    "reads.](figs/species-frac-B1242.pdf)\n",
    "\n",
    "[species-frac]\n",
    "\n",
    "The last question we want to answer is, “How many additional reads are\n",
    "needed to increase the fraction of observed species to a given\n",
    "threshold?”\n",
    "\n",
    "To answer this question, we need a version of `RunSimulation`\n",
    "that computes the fraction of observed species rather than the number of\n",
    "new species.\n",
    "\n",
    "```python\n",
    "# class Subject\n",
    "\n",
    "    def RunSimulation(self, num_reads):\n",
    "        m, seen = self.GetSeenSpecies()\n",
    "        n, observations = self.GenerateObservations(num_reads)\n",
    "\n",
    "        curve = []\n",
    "        for k, obs in enumerate(observations):\n",
    "            seen.add(obs)\n",
    "\n",
    "            frac_seen = len(seen) / float(n)\n",
    "            curve.append((k+1, frac_seen))\n",
    "\n",
    "        return curve\n",
    "```\n",
    "\n",
    "Next we loop through each curve and make a dictionary, `d`,\n",
    "that maps from the number of additional reads, `k`, to a list\n",
    "of `fracs`; that is, a list of values for the coverage\n",
    "achieved after `k` reads.\n",
    "\n",
    "```python\n",
    "\n",
    "    def MakeFracCdfs(self, curves):\n",
    "        d = {}\n",
    "        for curve in curves:\n",
    "            for k, frac in curve:\n",
    "                d.setdefault(k, []).append(frac)\n",
    "\n",
    "        cdfs = {}\n",
    "        for k, fracs in d.iteritems():\n",
    "            cdf = thinkbayes.MakeCdfFromList(fracs)\n",
    "            cdfs[k] = cdf\n",
    "\n",
    "        return cdfs\n",
    "```\n",
    "\n",
    "Then for each value of `k` we make a Cdf of\n",
    "`fracs`; this Cdf represents the distribution of coverage\n",
    "after `k` reads.\n",
    "\n",
    "Remember that the CDF tells you the probability of falling below a given\n",
    "threshold, so the *complementary* CDF tells you the\n",
    "probability of exceeding it. Figure [species-frac] shows complementary\n",
    "CDFs for a range of values of `k`.\n",
    "\n",
    "To read this figure, select the level of coverage you want to achieve\n",
    "along the $x$-axis. As an example, choose 90%.\n",
    "\n",
    "Now you can read up the chart to find the probability of achieving 90%\n",
    "coverage after `k` reads. For example, with 200 reads, you\n",
    "have about a 40% chance of getting 90% coverage. With 1000 reads, you\n",
    "have a 90% chance of getting 90% coverage.\n",
    "\n",
    "With that, we have answered the four questions that make up the unseen\n",
    "species problem. To validate the algorithms in this chapter with real\n",
    "data, I had to deal with a few more details. But this chapter is already\n",
    "too long, so I won’t discuss them here.\n",
    "\n",
    "You can read about the problems, and how I addressed them, at\n",
    "<http://allendowney.blogspot.com/2013/05/belly-button-biodiversity-end-game.html>.\n",
    "\n",
    "You can download the code in this chapter from\n",
    "<http://thinkbayes.com/species.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "The Unseen Species problem is an area of active research, and I believe\n",
    "the algorithm in this chapter is a novel contribution. So in fewer than\n",
    "200 pages we have made it from the basics of probability to the research\n",
    "frontier. I’m very happy about that.\n",
    "\n",
    "My goal for this book is to present three related ideas:\n",
    "\n",
    "-   **Bayesian thinking**: The foundation of Bayesian\n",
    "    analysis is the idea of using probability distributions to represent\n",
    "    uncertain beliefs, using data to update those distributions, and\n",
    "    using the results to make predictions and inform decisions.\n",
    "\n",
    "-   **A computational approach**: The premise of this book\n",
    "    is that it is easier to understand Bayesian analysis using\n",
    "    computation rather than math, and easier to implement Bayesian\n",
    "    methods with reusable building blocks that can be rearranged to\n",
    "    solve real-world problems quickly.\n",
    "\n",
    "-   **Iterative modeling**: Most real-world problems\n",
    "    involve modeling decisions and trade-offs between realism and\n",
    "    complexity. It is often impossible to know ahead of time what\n",
    "    factors should be included in the model and which can be abstracted\n",
    "    away. The best approach is to iterate, starting with simple models\n",
    "    and adding complexity gradually, using each model to validate the\n",
    "    others.\n",
    "\n",
    "These ideas are versatile and powerful; they are applicable to problems\n",
    "in every area of science and engineering, from simple examples to topics\n",
    "of current research.\n",
    "\n",
    "If you made it this far, you should be prepared to apply these tools to\n",
    "new problems relevant to your work. I hope you find them useful; let me\n",
    "know how it goes!\n",
    "\n",
    "[^1]: Based on an example from\n",
    "    <http://en.wikipedia.org/wiki/Bayes'_theorem> that is no longer\n",
    "    there.\n",
    "\n",
    "[^2]: Ruggles and Brodie, “An Empirical Approach to Economic\n",
    "    Intelligence in World War II,” *Journal of the American\n",
    "    Statistical Association*, Vol. 42, No. 237 (March 1947).\n",
    "\n",
    "[^3]: Zhang et al, Distribution of Renal Tumor Growth Rates Determined\n",
    "    by Using Serial Volumetric CT Measurements, January 2009\n",
    "    *Radiology*, 250, 137-144.\n",
    "\n",
    "[^4]: If you are not familiar with Python generators, see\n",
    "    <http://wiki.python.org/moin/Generators>."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
